// Copyright (c) Lawrence Livermore National Security, LLC and other VisIt
// Project developers.  See the top-level LICENSE file for dates and other
// details.  No copyright assignment is required to contribute to VisIt.

// ************************************************************************* //
//  File: avtToroidalPoloidalProjectionFilter.C
// ************************************************************************* //

#include <avtExtents.h>
#include <avtToroidalPoloidalProjectionFilter.h>
#include <avtVector.h>

#include <DebugStream.h>

#include <vtkCellData.h>
#include <vtkCell.h>
#include <vtkPointSet.h>
#include <vtkPoints.h>
#include <vtkUnsignedCharArray.h>

#include <limits>

// ****************************************************************************
//  Method: avtToroidalPoloidalProjectionFilter constructor
//
//  Programmer: hari -- generated by xml2avt
//  Creation:   Tue Feb 1 09:27:46 PDT 2011
//
// ****************************************************************************

avtToroidalPoloidalProjectionFilter::avtToroidalPoloidalProjectionFilter() :
  rMax(0)
{
}


// ****************************************************************************
//  Method: avtToroidalPoloidalProjectionFilter destructor
//
//  Programmer: hari -- generated by xml2avt
//  Creation:   Tue Feb 1 09:27:46 PDT 2011
//
//  Modifications:
//
// ****************************************************************************

avtToroidalPoloidalProjectionFilter::~avtToroidalPoloidalProjectionFilter()
{
}


// ****************************************************************************
//  Method:  avtToroidalPoloidalProjectionFilter::Create
//
//  Programmer: hari -- generated by xml2avt
//  Creation:   Tue Feb 1 09:27:46 PDT 2011
//
// ****************************************************************************

avtFilter *
avtToroidalPoloidalProjectionFilter::Create()
{
    return new avtToroidalPoloidalProjectionFilter();
}


// ****************************************************************************
//  Method:      avtToroidalPoloidalProjectionFilter::SetAtts
//
//  Purpose:
//      Sets the state of the filter based on the attribute object.
//
//  Arguments:
//      a        The attributes to use.
//
//  Programmer: hari -- generated by xml2avt
//  Creation:   Tue Feb 1 09:27:46 PDT 2011
//
// ****************************************************************************

void
avtToroidalPoloidalProjectionFilter::SetAtts(const AttributeGroup *a)
{
    atts = *(const ToroidalPoloidalProjection*)a;
}


// ****************************************************************************
//  Method: avtToroidalPoloidalProjectionFilter::Equivalent
//
//  Purpose:
//      Returns true if creating a new avtToroidalPoloidalProjectionFilter with the given
//      parameters would result in an equivalent avtToroidalPoloidalProjectionFilter.
//
//  Programmer: hari -- generated by xml2avt
//  Creation:   Tue Feb 1 09:27:46 PDT 2011
//
// ****************************************************************************

bool
avtToroidalPoloidalProjectionFilter::Equivalent(const AttributeGroup *a)
{
    return (atts == *(ToroidalPoloidalProjection*)a);
}


// ****************************************************************************
//  Method:      avtToroidalPoloidalProjectionFilter::GetCentroid
//
//  Purpose:
//      Gets the centroid of the phi=0 portoin of the mesh
//
//  Arguments:
//      inPts        Input point from the mesh
//
//  Programmer: hari -- generated by xml2avt
//  Creation:   Tue Feb 1 09:27:46 PDT 2011
//
// ****************************************************************************

avtVector avtToroidalPoloidalProjectionFilter::GetCentroid(vtkPoints* inPts)
{
    //get all the points, narrow down selection of points
    //that fall on plane at point 0,0,0 and vector 1,0,0 (yz plane)
    //and then pick points the that are on positive side of 0,1,0 (xz plane)
    //this should give me a circle of points..
    //add up the circle and compute the centroid and the distance from 
    //0,0,0 -> to this centroid is a R0
    //this is kind of hackish but should work for now..

    int nPts = inPts->GetNumberOfPoints();

    avtVector yzPlane(1,0,0); //yz plane
    avtVector xzPlane(0,1,0); //xz plane

    int count = 0;
    avtVector centroid(0,0,0);

    for(int i = 0; i < nPts; ++i)
    {
        //construct vector out of every point
        double pts[3];
        inPts->GetPoint(i,pts);

        avtVector pt(pts[0],pts[1],pts[2]);
      
        pt.normalize();

        //dot the two points, 0 means the point is on the XZ plane
        if(fabs(xzPlane.dot(pt)) < std::numeric_limits<double>::epsilon())
        {
            //std::cout << centroid << std::endl;

            //of these take only the positive side of the YZ plane 
            if(yzPlane.dot(pt) > 0)
            {
                //add up all these points and counts..
                //std::cout << centroid << std::endl;

                centroid += avtVector(pts[0],pts[1],pts[2]);
                count++;
            }
        }
    }

    centroid /= (double)count;

//    std::cout << "count: " << count << " " << centroid << std::endl;
    
    return centroid;
}
// ****************************************************************************
//  Method: avtToroidalPoloidalProjectionFilter::ExecuteData
//
//  Purpose:
//      Sends the specified input and output through the ToroidalPoloidalProjection filter.
//
//  Arguments:
//      in_dr      The input data representation.
//
//  Returns:       The output data representation.
//
//  Programmer: hari -- generated by xml2avt
//  Creation:   Tue Feb 1 09:27:46 PDT 2011
//
//  Modifications:
//    Kathleen Biagas, Tue Aug 21 16:28:11 MST 2012
//    Perserve coordinate type.
// 
//    Eric Brugger, Tue Aug 19 09:35:56 PDT 2014
//    Modified the class to work with avtDataRepresentation.
//
// ****************************************************************************

avtDataRepresentation *
avtToroidalPoloidalProjectionFilter::ExecuteData(avtDataRepresentation *in_dr)
{
    //
    // Get the VTK data set.
    //
    vtkDataSet *in_ds = in_dr->GetDataVTK();

    vtkDataSet *ds = in_ds->NewInstance();
    ds->ShallowCopy(in_ds);

    int do_type = ds->GetDataObjectType();

    if (do_type != VTK_STRUCTURED_GRID &&
        do_type != VTK_UNSTRUCTURED_GRID &&
        do_type != VTK_POLY_DATA)
       // silently fail for rectilinear grids
        return in_dr;

    vtkPointSet *ps = (vtkPointSet *) ds;
    vtkPoints *inPts = ps->GetPoints();

    int nPts = inPts->GetNumberOfPoints();

    vtkPoints *outPts = vtkPoints::New(inPts->GetDataType());
    outPts->SetNumberOfPoints( nPts );

//    double R0 = atts.GetR0();
//    double r = atts.GetR();

    double R0, r;

    avtVector centroid;

    if( atts.GetCentroidSource() == ToroidalPoloidalProjection::Manual )
    {
      centroid = avtVector( atts.GetCentroid() );
    }
    else //if( atts.GetCentroidSource() == ToroidalPoloidalProjection::Auto )
    {
      centroid = GetCentroid(inPts);
    }

    R0 = centroid.x;

    rMax = 0;

    double poloidal_angle_theta = 0.0, toroidal_angle_phi = 0.0;
    //test mapping of theta

//    double test_r = 1;
//    double angle = 45;

//    for(unsigned int i = 0; i < 360; i += angle )
//    {
//        //convert to radians
//        double test_theta = ((double)i) * (M_PI/180.0);
//        double test_u = test_r*cos(test_theta);
//        double test_z = test_r*sin(test_theta);

//        std::cout << "angle = " << i << " radian: " << test_theta << " z = "
//                << test_z << " u = "
//                << test_u << " atan(z,u)= "
//                << atan2(test_z,test_u)*180.0/M_PI << " atan(z,-u) "
//                << atan2(test_z,-test_u)*180.0/M_PI << std::endl;
//    }

    for (int i = 0 ; i < nPts ; i++)
    {
        double pt[3];

        //get point
        inPts->GetPoint(i, pt);

        //map to x,y,z to make it a little more readable
        double x=pt[0],y=pt[1],z=pt[2];

        //u = difference of major radius and projected radius of point
        //on xy plane
        double u = R0 - sqrt(x*x + y*y);

        //r = sqrt(u*u + (z-centroid.z)*(z-centroid.z));
        //compute hyp(r) = sqrt(length of adjacent(u)**2 + length of opposite(z)**2)
        //atan2 supports range of -PI to PI ..
        //theta = inverse tangent of opposite over adjacent
        //u is negated here because Linda Sugiyama requested the
        //mapping goes counter clockwise and outside of the torus
        //should be 0, inside is PI/-PI top is PI/2 and bottom is
        //-PI/2 and this mapping has that effect
        poloidal_angle_theta = atan2((z-centroid.z),-u);

        //phi is computed from and atan2 supports range of -PI to PI ..
        //x = (R0 +rcos(theta))sin(phi), y = (R0+rcos(theta)cos(phi))
        //x/sin(phi) = R0 + rcos(theta), y/cos(phi) = R0 + rcos(theta)
        //make equations equivalent -> x/sin(theta) = y/cos(theta)
        //tan(theta) = x/y -> theta = inverse tangent of (x/y)
        toroidal_angle_phi = atan2(x,y);

        double out_pt[3];

        //map to 2D x = toroidal, y = poloidal, z is projected to 0
        out_pt[0] = toroidal_angle_phi; //x axis goes in the direction of torus
        out_pt[1] = poloidal_angle_theta;//y axis goes in the direction of the poloid

        if( atts.GetProject2D() == true )
          out_pt[2] = 0.;
        else //if( atts.GetProject2D() == false )
        {
          r = sqrt(u*u + (z-centroid.z)*(z-centroid.z));

          out_pt[2] = r;

          if( rMax < r )
            rMax = r;
        }

        outPts->SetPoint(i, out_pt);
    }

    //set new points to dataset
    ps->SetPoints(outPts);
    outPts->Delete();

    //compute ghost cells to ensure bounds/wrap arounds are dealt with properly
    int ncells = ds->GetNumberOfCells();
    vtkUnsignedCharArray *gz = (vtkUnsignedCharArray *) 
                                 ds->GetCellData()->GetArray("avtGhostZones");
    if (gz == NULL)
    {
        gz = vtkUnsignedCharArray::New();
        gz->SetNumberOfTuples(ncells);
        gz->SetName("avtGhostZones");
        ds->GetCellData()->AddArray(gz);
        gz->Delete();
    }
    double b[6];
    for (int z = 0 ; z < ncells ; z++)
    {
         ds->GetCell(z)->GetBounds(b);

         //check bounds in x and y direction, if either x or y
         //direction have large bounds then presume that they are
         //wrapping around and therefore can be used as ghost cells
         if( (b[0] < -1 && 1 < b[1]) || (b[2] < -1 && 1 < b[3]) )
         {
             gz->SetValue(z, 1);
         }
         else
             gz->SetValue(z,0); //setting any other values to 0
    }

    // Update to get the correct rMax value.
    if( atts.GetProject2D() == false )
      UpdateDataObjectInfo();
    
    avtDataRepresentation *out_dr = new avtDataRepresentation(ds,
        in_dr->GetDomain(), in_dr->GetLabel());

    ds->Delete();

    return out_dr;
}


void
avtToroidalPoloidalProjectionFilter::UpdateDataObjectInfo(void)
{
    avtDataAttributes &outAtts     = GetOutput()->GetInfo().GetAttributes();

    outAtts.SetContainsGhostZones(AVT_CREATED_GHOSTS);

    if( atts.GetProject2D() == true )
      outAtts.SetSpatialDimension(2);

    else //if( atts.GetProject2D() == false )
      outAtts.SetSpatialDimension(3);

    GetOutput()->GetInfo().GetValidity().InvalidateSpatialMetaData();

    outAtts.GetOriginalSpatialExtents()->Clear();
    outAtts.GetDesiredSpatialExtents()->Clear();
    outAtts.GetActualSpatialExtents()->Clear();

    double bounds[6] = { -M_PI, M_PI, -M_PI, M_PI, 0, rMax };
    outAtts.GetThisProcsOriginalSpatialExtents()->Set(bounds);

    outAtts.SetXLabel("Toroidal");
    outAtts.SetYLabel("Poloidal");

    if( atts.GetProject2D() == false )
      outAtts.SetZLabel("R");
}
